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The extinction of a single species due to demographic stochasticity is analyzed. The discrete nature 
of the individual agents and the Poissonian noise related to the birth-death processes result in local 
extinction of a metastable population, as the system hits the absorbing state. The Fokker-Planck 
formulation of that problem fails to capture the statistics of large deviations from the metastable 
state, while approximations appropriate close to the absorbing state become, in general, invalid as 
the population becomes large. To connect these two regimes, a master equation based on a real 
space WKB method is presented, and is shown to yield an excellent approximation for the decay 
rate and the extreme events statistics all the way down to the absorbing state. The details of the 
underlying microscopic process, smeared out in a mean field treatment, are shown to be crucial for 
^ , an exact determination of the extinction exponent. This general scheme is shown to reproduce the 

, O ' known results in the field, to yield new corollaries and to fit quite precisely the numerical solutions. 

' Moreover it allows for systematic improvement via a series expansion where the small parameter is 

, , the inverse of the number of individuals in the metastable state. 
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Local extinction due to demographic stochasticity is a key issue in the analysis of persistence and viability of 
small populations. [U, [3, Q In particular, it allows ecologists to identify endangered species and to specify conservation 
policies, it dictates the appearance and disappearance of favored and neutral genetic mutations 4] and is of importance 
in the determination of the critical population size needed to support an epidemic. Q As in most cases of rare and 
O . extreme events the important quantities to be measured and compared with the theory are the expected extinction 
' time of the populations, [2] and the probability distribution close to the absorbing state. Technically, however, the 
' . effect of demographic stochastisity has traditionally been taken into account using some version of a Fokker-Planck 
I ^ , '■ equation, where the "diffusion coefficient" is a function of the population size.Q 

Recently, much interest has been focussed on the calculation of extinction rates for systems whose macroscopic 
dynamics exhibits a stable state which is nevertheless only metastable due to rare fluctuations which can drive the 
\ system to extinction. 0, H, Q In particular, it has been realized that in general the Fokker-Planck (FP) expansion 
about the (meta-)stable state is incapable of predicting the extinction rate. This is due to the fact that the Fokker- 
Planck expansion is only valid for up to 0{N'^^^) fluctuations to the large, 0{N), number of particles in the metastable 
state. The FP approximation fails to correctly describe the very large fluctuations necessary to reach the absorbing 
' state of zero particles. The FP treatment also smears out the microscopic differences between processes as it reflects 
\^ ' a local analysis close to the metastable fixed point. In order to get the correct statistics for rare and extreme events 
■ one should base the estimate on the exact Master equation that describes the stochastic process, and to employ the 
I method of extreme statistics, or more simply put, the WKB approximation, to solve the relevant master equation. 

Elgart and Kamenev Q made an interesting observation in this context: using the Peliti-Doi [Til, IT3| technique to 
map the exact master equation into a "quantum mechanical" problem (Schroedinger-like equation in imaginary time 
with second quantized Hamiltonian) they were able to identify the classical trajectory that connects the metastable 
' fixed point and the absorbing state. This identification allows them to calculate the classical ("geometrical optics") 
. ^ ' action along this trajectory, a first approximation to the extinction time. Assaf and MeersonQ then suggest a general 
spectral method to improve beyond the Elgart-Kamenev results, employing the generating function formalism and 
5h _ using the Sturm-Liouville theory of linear differential operators. 

In this paper, we will present a general scheme to deal with the local extinction problem, based on the time- 
independent "real space" WKB approximation (unlike Refs. ^, ^ who used a time-dependent momentum space 
presentation). The method presented is easy to use, its intuitive meaning is transparent, and its range of applicability 
covers, essentially, any single species problem. 

This paper will be organized as follows: In the next section we exemplify the technique for what is perhaps the 
archetypical problem in this class, a logistic birth-death process of a single species. Beside its importance, the solution 
of this example demands the use of all the components of the technique - a Fokker-Planck solution applicable close to 
the metastable fixed point, small n approximation close to the absorbing state and a WKB solution that encompasses 
the FP regime and connects to the small n region. This model, thus, serves also as a nice pedagogical introduction. 
The third section deals with a similar birth death process, but when the number of offsprings at each birth event is 
two, as in the case of domain walls in magnetic systems. Here, a series of mathematical "miracles" occur, which allow 
for a simple calculation of the extinction rate (for the case of an initial even number of particles) without recourse to 
the WKB method. In the fourth section the effect of a single agent death term is incorporated, and in the last section 
the marginal case of neutral mutation is analyzed. We then conclude with a summary and some final observations. 
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I. STOCHASTICITY, LOGISTIC GROWTH AND EXTINCTION 



In this section we study the stochastic dynamics of a combination of two fundamental processes: particles giving 
birth to new particles at rate a; and pair annihilation at rate (3. In that case the average number of particles is about 
a/ j3. The full "physical optics" solution for both the probability distribution, P„, of the metastable state and the 
extinction rate, in the limit where /3 ^ a is given, based on a WKB approximation for the exact master equation. 
We confirm our calculations by comparison to a direct numerical solution of the master equation. 

The microscopic rules that govern this process are: 

P ^ 2P 

P + P ^ (1) 
The exact master equation for P„, the probability of having n particles, is 

P„ = a [-nF„ + {n- l)P„_i] + | [-n(n - 1)P„ + [n + 2){n + l)P„+2] (2) 

At the mean-field level the process is described by the reaction equation, 

n — an — [3n{n — 1) (3) 

which has the stable solution n = a/ (3 + 1. Technically, this expression may be derived from the exact master 
equation ([2]) by calculating the time derivative of the average population (n) = ^ nP„ and using the approximation 
{n?) = (n)^. The stable solution becomes metastable due to the effect of stochastisity (in particular, all particles may 
annihilate each other and the system will be stuck in the absorbing state Pq = 1). Our aim is to calculate the typical 
time of this extinction event. 

Since that stochastic process has no memory (a Markov process) it may be described by a transition matrix that 
specify the rates to pass from one microscopic configuration to the other. Clearly, the absorbing state {Pq — 1, Pn — 
Vn 0} is an eigenvector of that matrix with an eigenvalue 0. All other eigenstates admit negative eigenvalues, and 
we denote the absolute value of the highest of these eigenvalues as F. The corresponding eigenvector is the stochastic 
metastable state, so our mission is to calculate F. Our main interest is in the case a/f3 ^ 1 so that the typical 
number of particles is large. This implies that the probability to reach the absorbing state is, as we will calculate, 
exponentially small. 

The metastability of the system implies that at long times the P„ decay exponentially as e~'"*. Thus we need to 
solve the master equation with the left-hand side replace by — rP„. However, since F is exponentially small, we can 
drop this term altogether. Technically this implies that we only have to solve for a steady state vector rather than 
doing time dependent semiclassical analysis which is much more complicated. 



A. Fokker-Planck equation and its limitations 



The standard approach for solving the now time-independent master equation is to transform it into a Fokker- 
Planck (FP) equation P, The nominal prescription for doing this is to expand P„±i, etc. in a Taylor series, 
dropping terms involving more than two derivative. However, phrasing the problem this way does not explain why, or 
more specifically, when this procedure is justified; moreover, the emerging equation is not unique - the same expansion 
may be done for nP„, for example, yielding different equation. The real justification underlying the Fokker-Planck 
approximation is that, for small P/a, Pn is a smooth function of the 0{1) variable 



y=\- 



1/2 



(3 



(4) 



Then, P„±i is equal to P{y± (3/q) and so may be expanded formally with regard to the small parameter 13 /a. The 
resulting series may be written as: 




CiP{y) + 



where 



d 



= 



(5) 



(6) 



3 



and 

£i = ^ (P'" + byP" + 8P' + AyP + 2y^P) (7) 

This perturbative expansion is justified (at least in the sense of asymptotic series) if, expanding P{y) = Po{y) + 
y/PaPi{y) + . . ., the correction term (3/aPi{y) <C Po{y)- Plugging this series into ([5]) and collecting terms order by 
order one finds that 

Then, looking at the equation for Pi 

CoPiiy) ^ CiPoiy), (9) 



one is able to identify that the leading correction is proportional to y^ y/ f3 / aPo{y) . Thus the the FP equation is only 
valid up to {n — a/ j3) ^ (a//3)^/'^, as mentioned above. This limit on the FP reliability is clearly demonstrated in 
Figure [U where for a metastable population of 100 individuals the FP solution is a good approximation to the exact 
probability distribution between 100 and 70. 

To find the decay rate, however, we need to have a solution valid down to n = 1, and the FP solution does not 
suffice. One can try to consider the low n limit of the master equation. 



B. Probability distribution close to extinction 

In the vicinity of the absorbing state one may use a simplified form of the master equation, exploiting the fact that 
P„ is a rapidly growing function of n. Then, P„ ^ Pn+2 and Pn-i ^ Pn, yielding the simplified recursion relation: 

^"+2 = f ^) , ^ " ^,, Pn (10) 

\ p J [n + 2)(n + 1) 

This approximate recursion relation leaves the even and odd n's decoupled, and has the solution 

_ /4ay-^ 2(fc-l)! 

with Pi, P2 arbitrary. Examining Eq. ^TU\i . we see that if Pi ^ P2, then as long as n <C a/f3, our assumption leading 
to the approximate recursion relation is valid. It is also clear that the rapid rise of the P„'s slows as n increases, with 
the Pn's reaching a maximum at n = 2a/(3. Clearly, however, the maximum probability state is at n « ck//?, so the 
recursion relation must fail before the bulk regime. In fact, our recursion relation works up to (3n/a ^ 1 whereas the 
FP solution only works when /3n/a « 1. There is no way to directly connect these two regimes. This is seen clearly 
in Fig. [U where the recursion relation and FP solutions are shown for the case P/a = 0.01. The resolution to this 
problem lies in the WKB method, which will allow us to connect the recursion relation results to the FP regime. 



C. WKB approximation and the extinction rate - the leading term 

To do WKB for our difference equation, [Tsl we write P„ = e'^", where Sn is assumed to be a smooth function of n, 
so that Sn±i ~ ± S!^. Since we already know (from the FP treatment) that the probability profile in the bulk is a 
Gaussian with width proportional to the square root of the metastable population, the quality of this approximation 
is controlled. Writing A„ = e'^", we have, assuming n ^ 1, 

= a (-n (-n^ + n^A^) (12) 
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FIG. 1: Probability distribution Pn/ P2 for (3 = 0.01, a = 1, for the basic model, Eq. ([TJ, together with the WKB approximation, 
Eq. (|30}, together with (fT6)) . (f29)) . and ((32]); the FP approximation, Eq. (f33|l : and the low/intermediate-n result, Eq. (fTTjl . 
Note that on the scale of the figure, the WKB approximation is indistinguishable from the exact result. 



or, simplifying, 

= /3nA^ - (2a + /3n)A„ + 2a 

= (A„-l)(/3nA2 +/3nA„-2a) (13) 

where we have factored out the trivial A„ = 1 root which is a result of conservation of probability. The other two 
roots for A are 



of which the larger, positive, one is relevant for us, since we want P„ to be an increasing function. This implies that 

(15) 

Integrating, we find that 




, , , 8a \ 1 / 8a 2a, f I3n (3n / 8^ , . ^ 1 

S-Q + nln + ^ - 1 + + S-ln 1- + 1 + + ^ -nln(2) + -n (16) 



(in J 2 \ Pn (3 \Aa 4a V (3n J ' ' 2 

The first important point to notice is that a/{(3n) extrapolates, in the interesting region, from 00 as n —> to unity, 
where n a/f3. In the first case 5' scales logarithmically with a/{(3n), hence Sn is proportional to a /{fin) and is 
large. In the second regime S" is almost constant and Sn scale with n which is also, in that case, large. Thus all the 
way to extinction Sn is large, as expected. Evaluating AS* = Sa/p — Sq, one obtains 

A5=^(l-ln(2)) (17) 
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in agreement with the result of Elgart and Kamenev. p 

We can make the connection to the formahsm presented in [6l| even more explicit, if we express the relations in 
terms of n(A) as opposed to A„. From Eq. (flSjl. 



thus 



AS* ^ I ln(A„)dn 





ln(A)^dA 
^ ^dA 

.ln(AHA)|r+^°"^dA 

-V^dA. 19) 
A 

If we now introduce the "momentum" p = 1/A and the "coordinate" q = nA, the expression for AS' may be rewritten 
as 

AS= [ qdp (20) 

where 







precisely reproducing Elgart-Kamenev equations for the action and the semiclassical escape path. The physical 
meaning of the "momentum" p is now clarified: it is the inverse of the geometrical growth rate of the quasistatic 
probability distribution. 

We now need to confirm that there exists an overlap region between the WKB solution and the n <^ a/ (3 recursion 
regime. For n a/ (3 our WKB solution for Sn may be approximated by, 

S„.So + in(l + ln(|^)) (22) 

To compare this with the recursion relation results, Eq. (jlip . in the limit n ^ 1: 



l^^peveu) ^ 1 J (^1 + In (^1^ ] ] + r ln(2) _ In ( ^ ) + ln(P2) 



2-v v.nyy 1'°<^>-'°(^) 

Indeed, the leading order asymptotics agrees. 

In t he ot her extreme our WKB result coincides with the FP treatment. Expanding the WKB solution for small 
y = yjl3/a{n - a/ (3) yields 

ln(S„) = ln(AS - 2/V3) (24) 

indeed reproducing the FP solution. 

D. WKB approximation and the extinction rate - first order corrections 

The explicit real space, time- independent WKB analysis allows one to go beyond the Elgart-Kamenev "geometrical 
optics" results and to obtain the leading corrections. In the previous subsection the generic substitution P„ = exp(S'„) 
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was implemented and is justified ex post facto by the fact that Sn tmns out to be 0{a/ (3) » 1. To proceed let us 
assume that 



P„ = cxp(5o(n) + Si{n) + S2{n) + ...) 



(25) 



where So{n) is the leading order Sn found above and Sm{n) is assumed to be 0{{f3 / a)'^~^). Beginning with the 
growth part of the master equation Tq, = a[— nP„ + (n — l)P„_i], plugging in Eq. l[25l) and expanding the small 
(©(/S/a)) terms in the exponent one has (note that any derivative adds a (/3/a) factor): 



ae 



So(n)+Si(«) + S2(n) 



5'o+5'i+S'2 




5q Si 



(26) 



While the first, 0(1), term in the bracket was used to determine A, the second 0{(3/a) term will be used here in 
order to find the function Si. Repeating that procedure and collecting the leading corrections from the annihilation 
part of the Master equation one finds. 



= 



nil 



Things simplify if we write S'l as a function of A: 



dA 



Si 



1 

2A 



(3n 
~2 



2/3nA^ 



2klS{ 



2n 



2KISI 



-A 



2a 
T 



(3n + 3f3nA^' 



(27) 



(28) 



so that all the coefficient functions are rational functions of A. The solution of this equation is best expressed in terms 



of Q„ 



pSi(n). 



/^(A„ + 1)^ 
V2A„ + 1 



so that to this order 



(29) 



(30) 



We see that Qn diverges as as n — > 0, since A„ diverges there, and vanishes as n~^/^ for large n, where A vanishes. 
Interestingly enough, there is no turning point, and the WKB approximation is good everywhere. This holds despite 
the fact that the coefficient of vanishes at A = 1, where S' vanishes, as is typical for WKB problems. In our case, 
the right hand side also vanishes at A = 1, so Q is regular there. We suspect that this is a consequence of the effective 
vanishing of F, so that the FP equation admits a trivial first integral, and so is effectively of first order. 

With this result for Qn in hand, we can finish the matching procedure. For n <C a/ (3, A„ is large as we noted, and 



A 



\/2a 
~i3n 



(31) 



Comparing to the recursion relation results, we have 



P2 



■2a2 



= Pi 



(32) 



This fixes the ratio of Pi to P2, which is precisely that which makes ln(P„) a smooth function. About the maximum, 
Kn~ 1, SO that Qn ~ 4A/%/3, and 



P„ 



V3 V" 



P,e^^e-^'/3 



(33) 
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FIG. 2: Decay rate vs. 1//3, for the basic model, Eq. ((T)), with a = 1, together with the analytic approximation, Eq. (|35|l . 
Here and in all other figures, the numerical calculation of the decay rate was performed using a quadruple-precision version of 
the sparse matrix eigenvalue solver ARPACK, applied to the master equation transition matrix, with the absorbing state 
eliminated. 



The sum over P„ is dominated by this FP Gaussian, and so, replacing the sum by and integral, we get 



(34) 



The decay rate F is then 



r - rJ-1- ~ , /i!!_p-2"(l-ln(2))//3 



(35) 



This result is plotted in Fig. [21 where we see the agreement is excellent for small (3. As F varies by so many orders of 
magnitude over the scale of the graph, it is impossible to see from this the role of the prefactor. In Fig. [3l we plot the 
ratio of F to e'^'^ as a function of f3. We see the accuracy is quite good and gets better as P gets smaller, as expected. 
This prefactor is, it should be noted, quite different from the factor a conjectured in Ref. [y]. In Fig. [H we present 
the graph of P„ together with the WKB, Fokker-Planck, and low/intermediate n approximations. We see that the 
WKB approximation is excellent everywhere, whereas the other approximations are more limited in their range of 
validity. In particular, the Fokker-Planck results is a serious overestimate of of the low n probability, and an equally 
bad underestimate of the large n probability. As the WKB and exact results cannot be distinguished on the scale of 
the figure, in Fig. |3]we present the ratio of the WKB to the exact result ioi a — 1, f3 = 0.01. We see that the WKB 
approximation is good to the expected few percent level, with it degrading slightly for very small n. Given that the 
WKB approximation is a large-n approximation, that it does as well as it does at low-n is a undeserved present, and 
a consequence of the remarkable accuracy of Stirling's formula down to n = 1. 

For completeness, we note that it is possible to derive a large /3 expansion of F as well. This is easily done by 
truncating the master equation matrix and computing its determinant as a power series in a/f3, then solving for F 
order by order. One finds 



F ; 



/3 



10 



38 



242 
135^ 



(36) 



It seems clear that this series has at best a finite radius of convergence. Both large and small f3 limits are compared 
to the exact results in Fig. \5\ 
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FIG. 3: Ratio of F to e^®, for the basic model, Eq. ([T}, compared to the WKB prediction, as a function of /3, for a — 1. 




FIG. 4: Ratio of the WKB approximation for P„ to the exact result for the basic model, Eq. ([1} for /3 — 0.01, a — 1. 



II. PARITY-PRESERVING MODEL 



In this section we will discuss the case where birth is to "twins", so that the even-odd parity of the number of 
particles is preserved. Here, a series of mathematical "miracles" occur, which allow for a simple calculation of the 
extinction rate (for the case of an initial even number of particles) without recourse to the WKB method. Furthermore, 
we present a solution of the probability distribution to all orders in (3 /a, so that the corrections are exponentially 
small. We also calculate the extinction rate to all orders in perturbation theory and manage to resum this divergent 
asymptotic series to obtain results correct to within exponentially small terms. 
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FIG. 5: Decay rate vs. 1/(5 for the basic model, Eq. ([T}, with a = 1, together with the large (3 approximation, Eq. (I36p . and 
the small /? approximation, Eq. p5|) . 



The model, 

P 3P 

P + P ^ (37) 

conserves the even/odd parity of the number of particles, so if the system is initiahzed with an odd number of particles 
it can never go extinct, instead reaching a steady-state. If the system is initialized with an even number of particles, on 
the other hand, the system can go extinct, with the survival probability decaying again as e"'"*, with T exponentially 
small for small /3. It should be noted that the mean-field equation for the model is the exact same as that of the 
original model above. 

We again start with the master equation, which now reads: 

Pn = -rP„ = I [-nPn + (n - 2)P„_2] + ^ [-n{72 - l)Pn + {n + 2){n + l)P„+2] (38) 

As above, T is exponentially small, and we may drop this term altogether. We again tackle the master equation by 
exploiting the fast growth of the P„'s for not too large n. This observation allows us to drop the second a term and 
the first (3 term, yielding the recursion relation 

\p J {n + l)(n -I- 2) 

which, up to a factor of 2, is the same as the approximate recursion relation we previously encountered. Now, however, 
the parity conservation implies that only the even terms are nonzero, with the odd terms being exactly decoupled. 
The recursion relation has the solution 

In principle, we should have to match this solution to the WKB solution, as we did in the nonparity case. The first 
miracle we encounter is that in fact the solution Eq. (PD|) is accurate throughout the Fokker- Planck region. To see 



this, note that for n = a/ (3 + y^a/(3, y ^ 0{1), the asymptotic expansion of P„ is 

P„« V2f^Ve"/2'5p2e-^'/4 (41) 
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It is straightforward to verify that this is the solution of the Fokker-Planck equation: 

= "(^(2/^) + 2^^j (42) 

It should be noted for the record that while Eq. is an accurate representation of P„ from n = 2 till past the 

peak, the Fokker-Planck Gaussian is again only valid in the peak region. 
Thus, we can use Eq. ([iP]) to calculate Ptot = J2k ^2fc- We find 

a\''^^ (fc-l)!2''- ^ 3 „ a 



For the moment, what is important is the leading order asymptotics of this, which can be calculated directly by 
applying Laplace's method to the sum. This is equivalent to integrating the Gaussian, and gives 

P,,,^V2^(^^y\'^/^Pp, (44) 
We are essentially done. The rate of probability flux out to the absorbing state is f3P2, which equals TPtot- Thus, 



r = ^ « ^ = , ^e-^'^^ (45) 



Note again that this is much smaller than the naive Fokker-Planck answer, which is proportional to e""/'*''. 

We can actually proceed to compute the corrections to this formula. One source of corrections is using the asymp- 
totics of The full asymptotics of 2-F2 for large argument are very beautiful: 



This is obviously a divergent series, and alternatively represents a resummation of the series. This result is easily 
proven by substituting it in the Hypergeometric Differential Equation. 

The second source of the corrections is the corrections to the P„'s due to the terms we dropped in the master 
equation. The structure here is also strikingly beautiful. If we denote our zeroth-order approximation of P„ by P°, 
we find 



P4 = pM 1 



The general trend is obvious: 



a 

Pe = Pg" ( 1 + - + 3 



P2k = P^, (1 + Xi(2m- 1)!! {^)] (48) 



Plugging this into the master equation shows that this is an exact solution (for F = 0, of course). Now, up to 
exponentially small corrections, Ptot is just multiplied by the correction factor: 

Ptot = PL(l+E(2m-l)!!f^ 



rn — 1 



'''' 
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FIG. 6: Extinction rate F and the relative error between the approximation, Eq. (|50p . and the exact value as a function of 1//3 
for a = 1. 



The final result for T is 



1/2/3 



2^2(1, 1; |, 2; ^) 



(50) 



Even though this answer is a resummation of the full asymptotic series, it is nevertheless not exact. The correction 
terms however are exponentially small, (relative to the exponentially small extinction rate) going like e-"/'^. This 
can be seen in the following graph, where we plot the relative error, comparing to an essentially exact numerical 
calculation (using extended precision arithmetic in Maple). Thus, for example, the error for a/P = 100 is six parts in 
lO^M 

For completeness, we also briefly write down the WKB solution. Firstly, 



n(A) 



so that 



Also, 



AS- 



A 



-dA 



a 



(51) 



(52) 



The solution for Q„ is the simple result 



5„ = - [ln(a//3) + 1 - ln(n)] 



(53) 



(54) 



All this can be seen to agree with our recursion relation solution. In fact, it implies that the recursion relation solution 
is valid everywhere, even past the FP regime. 

One can again calculate the large j3 limit of the decay rate. The first few terms of this series are: 



/3- 



a 



Ylo? 



24964a^ 



5 875/3 21875/32 58953125/33 95798828125/3^ 



(55) 



It would appear likely that this series is actually convergent. In any case, together with the small (3/a results above, 
they cover the entire range of parameters, as can be seen in Fig. [71 
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FIG. 7: Decay rate F for the parity conserving model, compared to the small /3 result, Eq. (I50|l . and the large /3 power-series 
expansion, Eq. (I55|l . The agreement between these analytical expressions and the exact results is so good in their respective 
spheres of validity, that the deviations are only visible quite far outside these areas. 



III. GENERAL NONPARITY MODEL 



Let us extend, now, our original, non-parity preserving, model to include a third process, the spontaneous death of 
particles at a rate 7 < a. Such a process appears naturally in many systems, from populations of animals (where 7 is 
the death rate of an individual) to the spread of a disease (where it correspond to a recovery of an infected agent, like 
in the SIR model [5j). We present the physical optics solution for this case also, again confirming it by comparison to 
the direct numerical solution. 

Adding the spontaneous decay of particles: 

P ^ (56) 

to our basic model, Eq. ([1]), changes both the mean field and the fluctuations. At the mean-field level this is equivalent 
to a simple change of the effective growth rate, a of aeff = a — 7, but this scaling is not true anymore if the fluctuations 
are taken into account. The master equation now reads: 

P„ = a [-nPn + in- l)F„_i] + ^ [-n{n - 1)P„ + {n + 2)(7i + l)P„+2] + 7 [-nPn + (n + l)P„+i] . (57) 
We start this time with the WKB solution. As before, the WKB ansatz yields an equation for A„, namely 

= [2a i~A„ + l)+f3n (-A„ + A^) + 27 (-A„ + A^)] 

-2a + /3nA„(A„-|-l) + 27A„] (58) 



A» - : 
2A„ 



with the solution 
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Now, A approaches the finite hmit, a/7 as n goes to 0, in contrast to the previous cases. Thus, we cannot solve the 
low-71 recursion relation by assuming that the P„ are increasing very rapidly. Rather, now P is irrelevant for low n, 
and we have to solve the f] — recursion. This is readily solved, for example by generating function techniques, and 
yields 



Pr.. 



7 



n(a — 7) 



^1 1 

7 



Pi 



(60) 



Indeed, for large n Pn grows geometrically, with the ratio a/7, agreeing with the small n WKB. 

We can now proceed with the remainder of the WKB procedure. As before, it is more convenient to work with 
n(A), given by 



n(A) 



a — 7A 
I3J A(A + 1) 



We see that in the mean- field regime, A « 1, the entire 7 dependence is through aeff- 
the situation is more complicated. 
Now, as before, S„ is given by 



(61) 



Away from this limit, however. 



Sn 



So 



n{A) ln(A) 



a/7 



«/7 



n(A 



dA 



So + n ln(A) + 



2(a-7A) 
/3A 



A 

2(a + 7) 

/3 



In 



A + 1 
A 



- In 



a + 7 



In particular. 



AS = 



2(a — 7) ^ 2(a + 7) a + 7 



/3 



2a 



(62) 



(63) 



We exhibit AS" as a function of 7 for fixed a in Fig. [S] We see that AS" decreases with increasing 7, vanishing 
quadratically as 7 approaches the threshold value of a. Also interesting is the dependence of AS* as a function of 7, 
for fixed aeff- This is show in Fig. [9l We see that as 7 increases, at fixed aeff, A5' decreases, leading to a faster 
decay rate due to increased fiuctuations. For large 7, in fact, AS vanishes as aeff /{2f3j). 
We now are in a position to continue to the calculation of Q„ = e^^^"'. The equation is 



2 

■7 



1 Q// 

2 ° 



nA„ 1 + S[ 



.25^ 



which simplifies to, using S" = l/(An'(A)) 

d 



■ - l3nA^ 
This has the solution 



7A 



dA 



Si 



1 
2A 



+ 2f3nA^ + J A 



n'{A) 
2n 



2a 



(3n + 3(3nA^ + 2jA 



(64) 



(65) 



Qr. 



A- 



Va3A„(A„ 



(a - 7A„)^a(2A„ + l)-7A2 



(66) 



where we have inserted the factor vc? so the definition of A reduces to that used in the 7 = case. 

The last step is to match to the low-n recursion relation solution. For Pn <C 1, the WKB solution reduces to 



2yl^a(a + 7) /a^ 
Pn V7^ 



(67) 
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FIG. 9: 2PAS/aeff as a function of 7, for fixed aeg. 
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We now need to calculate Ptot- Near the stable point, n = {a — 7)//3, 



P„ = AJ-^^e-^e-^i-^r (69) 
y cia — 7 a — 7 



which gives 



' na 



3 A 



Ptot - 4W— e 

y /3 a - 7 



AS 



;P,e^' (70) 



a + 7 (a — 7)^ 



The total probability flux out of the system is 7P1 + f3P2 ■ Since P2 is of order 1 relative to Pi , the (3 contribution to 
the flux is negligible, and so 



7P1 I a + ^ (a — 7) 



2 



Ptot V 2a 



-e 



-AS 



(71) 



This result is tested against the exact numerical answer in Fig. 1101 The results are again quite good, with the 
quality decreasing as 7 approaches a (for fixed P) as expected due to the decreasing equilibrium number of particles. 
The astute reader will note that our expression for F reduces in the 7 — > limit to our 7 = result above, even 
though some intermediate expressions (for example, the probability flux) do not correspond. We also note that the 
maximum value of A is 01/7, which approaches 1 as 7 approaches the threshold value of a. Thus, near threshold, 
the Fokker-Planck equation and the WKB treatment coincide. Of course, this solution still has to be matched to the 
low/intermediate-n solution to obtain the correct prefactor. 

It should also be noted that Doering, et al.,f7| investigated a class of models wherein all transitions are single- 
particle transitions, as opposed to the models investigated herein, which include a 2 particle annihilation process. 
This entire class of models can also be easily treated via our WKB method, and yields identical results to those of 
Doering, et al. Writing the birth term in the master equation by q;„ and the death term by 7„, the WKB solution can 
be written down for a very broad class of models where an = (l//3)a(/3n) and 7^ = (l//3)/?(n), and a, /3 are smooth 
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functions. In this case, the system admits a macroscopic metastable state with a large number of particles for small 
(3. In particular, we find 

A„ = ^ (72) 

SO that 

Sn - / ln(A„)dn (73) 







The Qn factor is: 



(74) 



which, upon matching to the low/intermediate-n result gives (up to an obvious typographical error) the result in 
Doering, et al. The advantage of the WKB method is that it generalizes to the multi-particle transition case. 

IV. THRESHOLD CASE AND THE LIFETIME OF A NEUTRAL MUTATION 

In that section let us consider the threshold case, a = 7. This case corresponds to the dynamics of a neutral 
mutation and has recently become the focus of extensive research, mainly in connection with HubbeU's unified neutral 
theory of biodiversity and biogeographv. |10i] Here, as we already saw for the near-threshold case, the P„'s are smooth 
and allow for a Fokker-Planck treatment. However, in this case the decay rate F is not exponentially small, and 
so cannot be ignored. First, let us consider what happens in the absence of /3. This problem was worked out by 
Pechenik and Levine. [isj They find that the mean number of particles is conserved and the variance grows linearly 
in time. Furthermore, the system exhibits a power-law (1/i) convergence to the empty state, and not an exponential 
dependence. This is indicative of a scale invariance in the problem. The /3 term serves to break this scale invariance, 
and gives a well-defined scale for the number of particles (for those replicas which still survive) . The original Fokker- 
Planck equation (ignoring boundary terms) was 

The /3 process introduces two new terms into the FP equation, arising from the Taylor expansion of 
{—n{n — l)P{n) + {n + 2){n + l)P{n + 2)) in the master equation. The first of these terms is a drift term, 
j3-^{'n?P). This is responsible for the term in the mean- field equation. The second term is a diffusion term. If 
we assume (3 small, then the additional diffusion induced by (3 can be ignored, and we are left with only the drift 
term. The new FP equation now reads 

We have assumed the time dependence is exponential, and are looking for the smallest eigenvalue A. We shall see 
that A vanishes in the /3 — s- limit, consistent with the power-law behavior found by P-L. 

It is useful to transform the FP equation to Shroedinger form. The first step in this process is to change variables 
to a; = ^Jn. This yields the equation 

_ aP" - (— + 2/3a;3)P' - 8(3x^P ^ ATP (77) 

X 

The next step is a similarity transformation to eliminate the first derivative 

P = x-^/^e-^^'/^'^Q (78) 

yielding the Schroedinger equation 



' Ax^ a I a 
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Fokker-Planck 




FIG. 11: Rescaled probability distribution (3 ^^^P„ as a function of P^^^n for the case a = f = 1, f3 = 0.01, together with the 
numerical solution of the Fokker-Planck equation. 

Clearly, in the absence of (3 there is no bound state. Rather there is a continuum that starts at zero. The potential 
with /3 has a single minimum, which is negative. Nevertheless, the ground state energy is positive, yielding a decay 
rate. The scaling of the decay rate is clear; by rescaling y = {jijaf^l^x^ (3 and a disappear from the equation, with 
the decay rate scaling as \/a(3. As a dverti sed, we verify the vanishing of A with (3. The presence of j3 has set the 
scale of (the surviving) n's, namely ^Ja/ (3, which is large for small (3. Lastly, the 1/x^ nature of the a potential is 
clearly a result of the scale-free nature of the (3 = problem. To get the prefactor multiplying the y/]3, we have to 
numerically solve the rescaled Shroedinger equation, (i.e., Eq. (j79p with a = (3 = 1) yielding the result 

r = l.lllVa^ (80) 

The resulting scaled Pn is shown in Fig. 1111 together with the rescaled exact numerical solution of the master equation 
for (3 = 0.01. We see that P„ is strongly peaked at the origin, corresponding to the zero particle mean- field solution. 

The scaling with (3 we obtained is the same as found by Doering, et al. However, as opposed to the other cases 
examined herein, the prefactor is different as now the mean time to extinction is not simply the inverse of the decay 
rate. This is due to the fact that all the eigenvalues of the master equation scale as whereas the other cases 
exhibited a single exponentially small eigenvalue. 

V. CONCLUSIONS 

We have solved for the fluctuation-induced extinction rates of various models exhibiting a macroscopic metastable 
state. Our primary methodology is use the WKB method for difference equations to directly solve the master equation. 
This WKB solution then has to be matched to the low-n P„'s, since the WKB method is a large n approximation. 
This technique is quite general and straightforward to implement, and produces quite accurate results as long as there 
are not too few particles in the metastable state. It reproduces the Doering, et al. results for the case of general 
one-particle transitions and generalizes to higher-order transitions. We have also shown the unique mathematical 
properties of the even/odd parity conserving model, where we are able to generate the full asymptotic expansion for 
the decay rate to all orders, and even to resum this divergent asymptotic series. 

One interesting point which arises from this analysis is the sensitivity of the decay rate to the exact form of the 
microscopic dynamics. This is apparent in the very different dominant e^'^ for the case of the parity and non-parity 
cases. This is also the case when one compares the logistic model with spontaneous decay studied in Section [Till to 
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the same model where the colhsion process P + P ^ is replaced by P + P ^ P at twice the rate. Even at the level 
of the Fokker-Planck dynamics valid near the metastable state, the widths of the distributions are different, with the 
variance (3a — j)/2(3 being replaced by a/ (3. The values of AS" in the two cases are very different, where Eq. ([63|) 
for the two-particle annihilation should be compared to 

AS^ In — = — ^ _ 1 In - 81 

Jo \1 + Pnj (3 13 

Thus, while the two expressions agree near threshold, where the Fokker-Planck description is sufficient, as 7 approaches 
0, AS* 2a/P{l — ln(2)) for the two-particle annihilation case, as we saw in Section 1, whereas AS* a/P, almost 
twice as large, for small 7 in the P + P ^ P case. Parenthetically, it should be noted that the small 7 limit is very 
singular in this latter case, as naively F seems to diverge as for small 7 ^ /3. In general, for small (3 and 7, the 

decay rate can be shown to be given by 

as so indeed vanishes as 7 — > 0, since there is no extinction in this case. In situations where the Fokker-Planck equation 
is valid all the way down to n = 0, as we saw was the case at threshold, at least the problem is parameterized by only 
two parameters, the center and width of the Gaussian. However, in general, the entire function n(A) is involved in the 
calculation of A^*. This should have important implications for the study of extinctions in the ecological community, 
for example, where reliable microscopic models are difficult if not impossible to obtain. 
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